A holistic finite difi"erence approach models 
linear dynamics consistently 

AJ Roberts* 
February 1, 2008 

Abstract 

I prove that a centre manifold approach to creating finite difference 
models will consistently model linear dynamics as the grid spacing be- 
comes small. Using such tools of dynamical systems theory gives new 
assurances about the quality of finite difference models under nonlin- 
ear and other perturbations on grids with finite spacing. For example, 
the advection-diffusion equation is found to be stably modelled for all 
advection speeds and all grid spacing. The theorems establish an ex- 
tremely good form for the artificial internal boundary conditions that 
need to be introduced to apply centre manifold theory. When numer- 
ically solving nonlinear partial differential equations, this approach 
can be used to derive systematically finite difference models which 
automatically have excellent characteristics. Their good performance 
for finite grid spacing implies that fewer grid points may be used and 
consequently there will be less difficulties with stiff rapidly decaying 
modes in continuum problems. 

Maths. Subj. Class: 37L65, 65M20, 37L10, 65P40, 37M99 
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1 Introduction 



Following the introduction of holistic finite differences in |T^, we would 
like to investigate numerical models for the dynamics of a field u{x, t) evolv- 
ing according to a nonlinear reaction-diffusion equation such as Ut = Uxx + 
f{u,Ux). This particular class includes Burgers' equation, / = —uu^, and 
autocatalytic reactions, such as / = u{l — u). However, before attacking 
such nonlinear problems, here we restrict attention to proving that the new 
methodology accurately models the djTiamics of quite general linear pde's. 

Modern dynamical systems theory has had to date very little impact 
on classical numerical approximations. Indeed, the very first sentence in 
Garcia-Archilla et al 



10 



says "Finite-element methods seem not to have 
benefited as much as spectral methods from some of the recent advances 
in the Dynamical Systems approach to partial differential equations". The 
concept of inertial manifolds has been developed to capture the long-term, 
low-dimensional dynamics of dissipative pde's [^. However, most efforts 
to construct approximations to an inertial manifold have been based upon 
the global nonlinear Galerkin method of Roberts Foias et al Q and 
Marion & Temam 



TBI 



al [0 and Foias 



This is so even for the variants explored by Jolly et 
In contrast, the approach proposed here is based purely 
upon the local dynamics on small elements while maintaining, as do inertial 
manifolds, fidelity with the solutions of the original pde. 

1 propose to use centre manifold theory to construct finite difference 
models. For problems in one spatial dimension consider implementing the 
method of lines by discretising in x and integrating in time as a set of ordinary 
differential equations, sometimes called a semi-discrete approximation 0, |^, 
e.g.]. We only address spatial discretisation and treat the resulting set of 
ordinary differential equations as a continuous time dynamical system. Clas- 
sical finite difference approximations are made by appealing to consistency 
in the limit as the grid spacing /i — 0; traditionally one constructs models to 
errors O [h"^) or O [h^) depending upon small h asymptotics, shown schemat- 
ically by the rightward- arrows in Figure |l|. In contrast, we here analyse the 
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dynamics at fixed grid spacing h and use centre manifold tlieory to accu- 
rately model the nonlinear dynamics — tlieory p|, e.g.] assures us that the 
low-dimensional, numerical model then accurately captures the dynamics in 
an expansion in the nonlinearity, shown schematically as the forward- arrows 
in Figure |l|. The analysis rests upon the exponential decay of the small, 
subgrid structures in each local element. Being essentially local in space, the 
analysis here is flexible enough to subsequently cater for spatial boundaries 
and spatially varying coefficients. I call the model "holistic" because the cen- 
tre manifold is made up of actual solutions of the pde and is thus invariant 
under algebraic rewriting of the governing equations. However, to apply the 
centre manifold theory we have to use a homotopy in a parameter 7: when 
7 = the discrete finite elements of the domain are completely uncoupled 
from each other; when 7 = 1, the requisite continuity is reclaimed and the 
physical pde solved. The caveat is that the centre manifold model has to 
be used at 7 = 1 whereas the supporting theory only guarantees accuracy 
in a neighbourhood of 7 = 0; we aim to make the useful neighbourhood big 
enough to include the relevant 7 = 1 (this sort of technique has proven effec- 
tive in thin fluid flows 0, e.g.]). One way to reasonably secure the centre 
manifold model, and the way explored herein, is to require that the model 
is also consistent with the pde as the grid spacing h ^ 0. Thus we aim 
to construct finite difference numerical models that not only are justified by 
their asymptotic expansions in nonlinearity and 7, but are also justified by 
an asymptotic expansion in h (see Figure |ip. This dual justification is the 
completely novel feature of the approach. 

The first step, taken in this paper, is to establish a centre manifold ap- 
proach that is also guaranteed to construct a consistent finite difference model 
for a general linear pde, shown schematically in Figure |l| as the disc in the 
7/i-plane (zero nonlinearity). We leave to later research the problem of guar- 
anteeing the consistent modelling of nonlinear dynamics. Herein we explore 
the finite difference modelling of the linear pde 



where the linear operator A, presumed generally dissipative, is even (it con- 
tains only even order derivatives in x with constant coefficients), and B is 
an odd linear operator (it contains only odd order derivatives with constant 
coefficients); A is assumed dissipative for all modes except u = const. The 
case of space-time varying coefficients to the linear problem is also left for 
later study; however, expect that because the analysis here is local in x, then 
such varying coefficients can be treated as a perturbing influence to the basic 
analysis herein. As a specific example, we discuss in §^ the linear advection- 



— = Au + eBu , 
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nonlinearity 



c onsistency » phy si ral problem 



7=1 




consistency ^^^^^^^ ^-^^^^ 

problem 



Figure 1: conceptual diagram showing: the traditional finite difference mod- 
elling approaches (rightward-arrows) the physical problem (upper disc) via 
asymptotic consistency as the grid size /i — » (left circles); whereas the holis- 
tic method approaches (forward-arrows) the physical problem via asymp- 
totics in nonlinearity and the inter-element coupling 7 (from right circle). 
Herein we estabhsh how to use the hohstic approach to do both in order to 
model a general linear problem (lower disc). 
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Figure 2: schematic picture of the equi-spaced grid, Xj spacing h, the un- 
knowns Uj, the artificial internal boundaries between each element (vertical 
lines), and in the neighbourhood of Xj the field Vj{x, t) which extends outside 
the element and, if 7 = 1, passes through the neighbouring grid values Uj±i. 



diffusion equation ut = —eux + Uxx and discover many remarkable properties 
of the holistic finite difference models. The separation between the two types 
of linear terms into Au and eBu occurs because first we prove consistency for 
even terms, §§, before moving on to prove consistency for the odd terms, §|^. 
The results are verified for the perturbed diffusion equation by the computer 
algebra program listed in the Appendix. 

Introduce a regular grid as shown in Figure |^ with grid points a distance 
h apart, Xj = jh for example, and using Uj to denote the value of the field 
at each grid point, 

Uj{t) =u{xj,t) . (2) 

We express the field in the neighbourhood of the jth grid point hj u = 
Vj{x, t). We do not restrict the function Vj to just the jth element, but allow 
it to extend analytically out to at least the adjacent grid points as shown in 
Figure 0. 

Herein we establish the small h consistency that follows from using the 
nonlocal, internal "boundary conditions" 

Vj{xj±i,t) = (1 - l)vj{xj, t) + 7^i±i(a;j±i, t) ■ (3) 

That is, the field of the jth element when evaluated at the surrounding 
gridpoints, Vj{xj±i,t), is a continuation between two critical extremes: when 
7 = 1, it is the field at those grid points, Vj±i{xj±i, t), to in effect recover the 
physical continuity as shown schematically in Figure H; but when 7 = 0, the 
field is just identical to the mid-element value Vj{xj,t) so that the element 
becomes isolated from all neighbours. Equivalently, @ is transformed to the 
following appealing two difference equations^ evaluated at the centre of the 

^Natural symmetry also makes the general results most easy to express in terms of 
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element, x = xj, 



fJ'xSxVjix,t) = 'yfi5vj{xj,t) and 6^Vj{x,t) = 'yd Vj{xj,t) . (4) 

That is, in the two extremes: when 7 = the first and second differences 
have to be zero; whereas when 7 = 1 the first two differences of the field 
centred on each element have to agree with the first two differences of the 
grid values. Note a distinction which is very important throughout this 
work: unadorned difference operators, such as the central mean /i = 
i?~^/^)/2 and central difference S = E^/"^ — E'^^"^ written in terms of the shift 
operator Euj = Uj+i [0, p64,e.g.], apply to the grid index j (with step 1) 
whereas those with subscript x, as in /ij. and 6^, are differences in x only (with 
step h). Using the definition of the amplitudes (|^), these internal boundary 
conditions (iBC) simplify to the following form which we use throughout 
evaluated at a; = x.- 



fJ'xSxVj{x,t) = 'jfj.5uj and 5lvj{x,t) = 'y5'^Uj . (5) 
In actually developing finite difference models these iBC's may take any of 



many equivalent forms |21, e.g.]. Small h consistency seems easiest to estab- 



lish in this particular discrete form. 



2 The dynamics collapses onto a centre man- 
ifold 

I establish here the basis of a centre manifold analysis of the linear pde (p. 
It appears necessary, and is the route taken here, to separate the linear effects 
into those generated by even terms, represented by A and presumed generally 
dissipative on the grid scale h but with one eigenvalue corresponding to 
u = const, I and those generated by odd terms, represented by eB. The 
analysis is to be based on the situation when the coefficient of the odd terms 
e = and when adjacent elements are decoupled, 7 = (the right-hand circle 
in Figure p. Then centre manifold theory 0, e.g.] guarantees the existence 
and relevance of a numerical model parametrised by the discrete values of 
the field at the grid points, uj. The constructed model is accurate to the 
order of the residuals of the differential equation (|l|). 

centred difference operators and so I use them throughout this paper. 

^The analysis presented here could be applied to modelling unstable dynamics, such 
as that from negative diffusion ut = —Uxx- The difference is that the relevance theorem 
would no longer apply — Ai would not be attractive and the finite difference model would 
not capture the long term dynamics. The centre manifold model would, however, capture 
all the finite solutions, if any. 
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The spectrum of the hnear dynamics is used to show there exists a centre 
manifold. Adjoin to the trivial dynamical equations 

7 = e = 0, (6) 

so that terms multiplied by e or 7 become "nonlinear terms" in the asymp- 
totic expansion we develop about 7 = e = 0. Setting e = eliminates the 
odd terms to leave simply the linear equation Ut = Au; for example, the 
diffusion equation Ut = u^x- This is to be solved in the vicinity of Xj for 
a field u = e'^^w{x) where here the dependence upon j is implicit in the 
eigenfunction w of Aw = Xw. This is a constant coefficient ODE and so has 
trigonometric general solutions w oc e*"^ with corresponding eigenvalue A(a) 
which is negative for non-zero wavenumber a as ^ is presumed generally 
dissipative; for example, A = — for the diffusion equation. The appropri- 
ate boundary conditions come from the nonlocal decoupling conditions that 
w{xj±i) = w{xj) from (§) with 7 = 0. Thus within each element the eigen- 
modes are: exp[ia„(x — xj)] for even integer n, where a„ = nn/h; and also 
sin[a„(a; — Xj)] for odd integer n. The spectrum is then A„ = X{n7i/h); for 
example, A„ = — n^vr^/Zi^ for the diffusion equation. Thus linearly and in the 
absence of inter-element coupling, generally expect all modes to decay expo- 
nentially quickly to zero on a time scale O {h^) as A(a) will be symmetric, 
except for the neutral mode, n = 0, which is constant in x, Vj{x,t) = uj. 
Thus for small enough e and 7, theory ([Q, p281] or p96]) assures us that 
there exists a centre manifold Ai for the system (|l|) coupled across elements 
by (|^). The centre manifold Ai is here, by (|^), to be parametrised by the 
values of the field at the grid points, Uj. Thus using u to denote the set of 
grid values Uj, the "amplitudes", theory p|, ||, ^ supports our description 



of the centre manifold and the evolution thereon as 



Uj = v{u, x) , such that iij = g{u) , (7) 

where the fact that the right-hand side functions depend upon the element j 
is implicit (no confusion need arise because the translational invariance in x 
leads to identical expressions in each element except for appropriate changes 
of the subscript j). The evolution iij = g{u) forms the holistic finite dif- 
ference model. When the model is constructed to errors O (^7^^ , then we 
account for interactions among i — 1 elements on either side of any given ele- 
ment and so the resulting finite difference model has a stencil of width 2i — l 
on the spatial grid. I call these models "holistic" because, unlike traditional 
finite difference modelling which just analyses separately each term in the 
equations, here Ai is made up of actual solutions of the pde and so here the 
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discretisation models all the possible interactions between all the terms in 



the equations [0, 0, e.g.]. 

Moreover, the evolution on the centre manifold forms an accurate low- 
dimensional model of the dynamics of the pde (|1|). Again provided e and 7 
are small enough, theory, [|, p282] or ||2^, pl28], assures us that all solutions 
sufficiently near the centre manifold Ai are not only attracted to Ai but 
exponentially quickly approach the actual solutions of the pde that make 
up Ai. The rate of attraction is approximated for practical purposes by 
the leading negative eigenvalue; for example, Ai = — 7r^//i^ for the diffusion 
equation. In the development of inertial manifolds by Temam and others, 
this property is sometimes termed the asymptotic completeness of the model, 
for example see Robinson [^] or Constantin et al 0, §12-3], and sometimes 
as exponential tracking p, e.g.]. Observe that this is one of the crucial new 
aspects brought to finite difference modelling by centre manifold theory. It 
asserts that on a finite grid spacing we will faithfully track the solutions of 
the original pde. There will be some limitations: steep gradients and large 
nonlinearities will test the model as always. But the centre manifold theory, 
as seen here in §^ and in introductory work [^T], [l^, provides a rationale and 
a method to construct the requisite adjustments to a finite difference model 
to robustly model a wide range of dynamics on a finite grid spacing. The 
main limitation is the rate of attraction to M.: the centre manifold model 
should be able to capture any dynamics occuring on a time-scale larger than 
l/|Ai| [Y? j-K^ for diffusion). Thus the model evolution on A^, (|^), captures 
all the long term dynamics with some provisos. 



3 Advection-difFusion is modelled robustly 

Here we explore perhaps the simplest nontrivial example in the class of pde's: 
the advection-diffusion equation 

ut = -eu^ + u^^ , (8) 

where e is the advection speed that will be treated as small in the asymptotics 
but investigated over a range of sizes. This pde fits into the scheme outlined 
in the previous sections with the operators A = and B = —d^. We show 
consistency for small grid spacing and find interesting and stable upwind 
approximations for large eh. The associated sophisticated dependence upon 
e is perhaps indicative of the need to treat odd operators differently to even. 

The centre manifold models discussed here were all derived by the com- 
puter algebra program given in Appendix 0. The listing in the appendix 
replaces the recording of tedious intermediate algebraic steps. Suffice to say 
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that the algorithm is based upon an iterative refinement technique reported 
T9| : given an approximation to the centre manifold and the evolution 



m 



thereon (|^), the computer algebra calculates a correction to reduce the resid- 
uals of the governing equations, the pde the iBC's (§) and the amplitude 
condition @. Iteration then drives the residuals to zero to some order of 
error. Thus by the Approximation theorem & O, e.g.] we are sure that the 



model is correct to the same order of accuracy. In this section we proceed to 
use the results without any further description of the algebra involved. 

First explore the holistic finite difference model with only first-order in- 
teractions between adjacent elements, that is, the case i = 2 where errors are 
O (7^). As mentioned earlier, we find odd derivatives in x cause a hierarchy 
of refinements parametrised by increasing powers of e. To low-order in e the 
computer algebra shows the evolution on the centre manifold is 
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u. 



+o(e^7^) 



(9) 



substitute 7 = 1 to obtain a model for the advection-diffusion pde (|]). 
Introduced automatically in this analysis is the novel term involving the 
square of the advection speed One alternative is to view this term 

as an upwind correction to the finite difference approximation of the first 
derivative: 
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which increases the weight of the upwind grid point is upwind if e 

is positive). Such upwind corrections are well known to be stabilising for 
finite advection speeds e. The other alternative is to view the novel term as 
increasing the dissipation operator, to 
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1 + 
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and thus also improving the stability of the finite difference scheme for finite 
speeds e. It is this latter view that seems easiest to establish at higher order. 

It is interesting to explore in some detail to high order in e. I find the 
model is (after setting 7 = 1) 



fi6 5^ 



where ui 



1 + 



h ' 



(10) 
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This is equivalent to 



^2 

ut = -eu^ + Ur,:r + — (e - dxfua^x + O {fi^^ , (11) 

and so is indeed consistent to O {h?) as — > , independent of e, with the 
advection-diffusion pde (|^). The coefficient of the enhanced diffusion is the 
famihar value 1 when eh is small. However, when the advection speed is large 
enough (taking e > hereafter in this section for simplicity), the diffusion 
coefficient vi increases to aid stabilisation. The series for vi in (0) is identical 
to that for (e/i/2) coth (e/;./2), and numerical summation of the series using 
the Shanks transform JI], §8.1] suggests strongly that z/i ~ e/;,/2 as e/i — > cxd 
and is indeed within a few percent of this value for e/i > 4 , see Figure |^. 
Thus for large advection speed e on a finite width grid, the centre manifold 
analysis promotes the model0 (written in terms of the backward difference 
operator V) 

iij ——Vuj = — — (uj — Uj^i) . (12) 
h h 

This is not, and need not be, consistent with the pde (|^) as /i — because it 
applies for finite eh. That it should be relevant to (^ comes from the centre 
manifold expansion in 7 albeit evaluated at 7 = 1 (via the "holistic" arrows in 
Figure |l]). Exact solutions of ([T^) are readily obtained. For example, consider 
a point release from j = at t = 0: Uj{0) = 1 if j = and is otherwise. 
Then the moment generating function G{z, t) = Y^JLo ^'''^ji't) is easily seen to 
be that for a Poisson probability distribution with parameter et/h, namely 
G{z,t) = exp[(2; — l)et/h]. Hence the mean location and variance of uj are 

IJ-j = cr'i = -r =^ f^x = and = eht . (13) 
h 

Thus for eh not small: this model has precisely the correct advection speed e; 
and although the variance is quantitatively wrong, eht instead of 2t, at least 
it is qualitatively correct for finite eh. The centre manifold model (|TD|) is 
consistent for small h and has the virtue of being always stable and will always 
maintain non-negativity of solutions no matter how large the advection speed 
e. 

Second, explore the holistic finite difference model with second-order in- 
teractions between adjacent elements, that is, the case i = 3 for which errors 



■^If the advection velocity is negative e < 0, then various signs change and the large 
eh model remains an upwind model, but is consequently written in terms of forward 



differences. This also occurs for the later model (nju accurate to errors O (7 
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Figure 3: coefficients of the centre manifold models ([ToD and ([T5| ) as a func- 
tion of advection speed and grid spacing, eh. These curves are at least of 
graphical accuracy and are obtained via the Shanks transform of the Taylor 
series to errors O (e 
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are O (7^), but to high order in the advection e. Computer algebra derives 
the model 
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h 

+ 0(7^e^°)• (14) 

Observe the marvellous feature that when we evaluate this at 7 = 1 all the 
terms in 6"^ associated with high orders of eh neatly cancel. The model for 
the advection-diffusion thus reduces to 
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See that in this model, as h 



2520 75600 2395008 
0, the hyperdiffusion coefficient 1^2 ~ 1/12 and 
the dispersion coefficient ^2 ~ 1/6 to give the classic second-order in h cor- 
rections to the central difference approximations of the first two derivatives. 
Indeed the model (ITl) is equivalent to 



Ut 



-(e - a 

90^ 



K O [K') , (16) 
, independent of e, with the advection- 



and so is consistent to O (/i^) as h 
diffusion pde (g). 

For large advection speed or grid size, large e/i, the model (|15D is aston- 
ishingly good. Using the large eh approximations indicated by the Shanks 
transforms plotted in Figure |^ for z/2 and K2, the model (^) reduces to simply 
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{Uj-2 - 4:Uj-i + 3Uj) + — {Uj-2 
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(17) 
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This large eh model uses only backward differences to incorporate second- 
order upwind estimates of the derivative, V + |V^, and the second derivative, 
V^. To show its good properties,^ consider again a point release from j = 
at time t = 0. The moment generating function G{z,t) 
the evolution governed by (0) is readily shown to be 

et , t , 



Y.T=oZ'u,{t) for 



G{z, t) = exp 



Then since 
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dz 
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and 
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we determine the mean position and variance of the spread in Uj to be 



et , 2 2t 
N = J and a,- = — 



Ail 



et and = 2t . 



(19) 



This predicted mean and variance following a point release are exactly correct 
for all time for the advection-diffusion pde (H). This specific result applies 
to all finite advection speeds e and all finite grid spacings h whenever eh is 
large enough. 

It seems that creating finite differences which, as shown in Figure 0, are 
both consistent for small grid spacing h and also holistically derived via centre 
manifold theory thus can lead to models which are remarkably accurate and 
stable over a wide range of parameters. The hierarchy of refinements in the 
coefficient, e, of odd derivatives, here just dx, seem to be useful for robust 
performance for finite h. 



4 Models of the even terms are consistent 

In this section I prove that the proposed centre manifold approach consis- 
tently models all the even derivatives in A. That is, the equivalent pde of 
the finite difference model on the centre manifold M. matches Ut = Au to 
some order in grid spacing h; the order of error is controlled by the order 
of truncation, i, in the coupling coefficient 7. The veracity of the following 
theorems are supported by the results of suitable variants of the computer 
algebra program of Appendix ^. 

Remarkably, the polynomials found here to describe the structure within 
each element are independent of the linear operator A\ 

^The upwind difference model ( pT| ) is only stable for e/i > 2/3 . However, from Figure |3| 
the approximation ( p^ ) is only applicable to ( p^ ) for eh greater than about 4; thus its 
instability for very much smaller eh is irrelevant. 
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Theorem 1 The centre manifold model 1^ constructed with the amplitude 
condition the internal boundary condition (Qj and to errors O (^7^ j , 
forms a sem-discrete finite difference approximation to the pde Ut = Au 
consistent to O (<9^^), where A is any even operator. 



Proof: I construct the proof in stages beginning with the end result and 
finding successive sufficient conditions for the preceding steps. Observe that 
I actually prove a slightly stronger result: by allowing the even operator A to 
contain a constant term oq, see (p7|), I prove the consistency of an invariant 
manifold model based upon the mode with eigenvalue Aq = cio- When ao > 
this forms a centre-unstable manifold model for which a relevance theorem 
also ensures asymptotic completeness. 

• To solve Ut = Au to errors O {^^^ , expand the centre manifold model 

to O (7^) as 

i-i £-1 
u = v{u, x) = Y^ 7^f ''(n, x) , and uj = g{u) = ^ ^^g'^Uj , (20) 

fc=0 fc=0 

where are some difference operators, and where v and g implicitly 
refer to the jth element. Note that superscripts to 7 denote expo- 
nentiation whereas those on v and g denote the index of coefficients 
in the asymptotic expansion: I do this because v and g have an im- 
plicit subscript j denoting the grid element under consideration. Then 
substitution into the pde and extracting powers of 7 shows that we 
require 

Av^ = g^v^ + g^v^-^ + ■■■ + g^-^v^ + g^v^ for < A; < ^ . (21) 

Similarly, substitution into the amplitude condition @) requires 

v^{u,Xj) = Uj, and v'^{u,Xj) = Q ioi I < k < I . (22) 

Whereas substitution into the IBC (^ requires, evaluating the left-hand 
sides at Xj, 

,A„' = {f-*;;; (23) 

and = {f;- (24) 

Equations (PT|-P^ form a well-posed system of equations for and g^ . 
In many applications of centre manifold theory, because A — g^ is sin- 
gular, we often solve each level in the hierarchy of equations in two 
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steps: the first is to find by ensuring that the right-hand side of ^ 
is in the range of ^ — , this is the so-called "solvability condition"; 
the second step is to find Vk- However, here we proceed to construct 
straightforwardly the solution of the entire set of equations in general. 



We show the hierarchy of differential equations (|2l|) are satisfied by 



functions and give consistent finite difference operators g^, if 
satisfy the recursive difference equation 

5lv^ = 5'^v^~^ for all x, and = constant . (25) 

— By the following induction argument (EH) implies that 



s:2m k I 5'^"^v'' ™ , for m = 0, . . . , , 



^ 0, m = k + l,k + 2, . . . . ^ ' 



Now, 



^2^2(m-l)^A:-l ^ COmmutC 

^2g2{m-l)^k-m ^ ^ ^^^^^ for m - 1 

^2m^k—m 



Then since (p6|) is trivially true for m = 0, it follows by induction 
that (p6D holds for all m provided m < k. Further, since is 
constant by (|25|), 6l''v^ is constant, so higher order differences 
(m > k) are all zero. 

By an "even" operator A 1 mean one which only involves even 
order derivatives in x. Hence write A formally as an infinite sum 
of even powers of the central difference operator 

oo 

a=y: a-^r , (27) 

m=0 

for some coefficients am', for example, for the diffusion operator 
in d), from [0, p65,e.g.], 

^ = ± arcsinh^d^.) = - ^^5^ + - . . . ; 

more generally, A could be any symmetric convolution operator 
for which the infinite sum (E^ forms a reasonable representation. 
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Then (pG]) ensures (^T]) since 

oo 
m=0 



m=0 

5f%^ + g^v^^'^ H h g''~^v^ + 5f V 



provided (7^^ = atS"^^ which are precisely the operators required to 
make the model Uj = g{u) of Ut = Au consistent to O (<9^^), when 

truncated as in to errors O (^7^^ . 

By simple substitution, a sequence of functions satisfying the re- 
currence (|25D , amphtude conditions ( ^2]) and the internal boundary 
conditions (53-121) are 



v'' = uj, v"" =pk{i)li5^^-\ + qk{i)5^^Uj, forfc>l, (28) 

where, as always, ^ = {x—Xj)/h is a grid scaled coordinate and provided 
that for k > 1 

6lpk=Pk-i, Pk{±k) = ±l, and pfe(0 = Ofore = 0,±l, (29) 
and similarly 

^iQk = qk-i , qk{-^k) = +1 , and qk{0 = for ^ = 0, ±1 , (30) 
after defining go (a;) = 1 and po{x) = 0. 

Analysing the difference tables for p^ and q^ and straightforward in- 
duction using (p^-|30|) proves that Pk{0 = Qk{0 = ^'^^ integer ^ G 



[— A; + 1, A; — 1]. Then the following pk and qk are the unique polynomi- 
als, of degree 2A; — 1 and 2k respectively, also satisfying pk{±k) = ±1 
and qki±k) = +|: 

1 '"'^ ^ 

PkiO = ^2k-l)\ n - ' and ^^(0 = ^PkiO , (31) 

as plotted in Figure For example, Pi(0 = ^ and ^1(0 = 

These polynomial pk and qk also need to satisfy the recurrences in (PU|- 



30) pointwise in ^. This is trivially true for k = 1. Now, 5^pk+i is 



from (^) a polynomial of degree 2k — 1, from its difference table has 
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Figure 4: graphs of the polynomials ( pT]) forming a basis for the fields of the 
approximations to the centre manifold. 
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the same zeros as pk, it is ±1 at ±k, and so must be Pk{0 
Similarly for 

4k 

Corollary 2 It follows immediately from the theorem that if the highest 
derivative in the even operator A is of order n, then the finite difference 
model is consistent with Ut = Au to O {h'^'''~^^ as h ^ 0. 

5 Odd perturbations are also consistent 

The results of the previous section on the modelling of even operators A are 
extremely satisfactory. The modelling of odd operators is not quite so neat. 
In the general linear pde (|1|) I introduced the odd terms, B, with a multiply- 
ing e. The reason is, as seen in that such odd terms generate extra terms 
in the finite difference model which are nonlinear in the coefficients of the odd 
derivatives, that is 0{e^). As ellaborated in §^ for the advection-diffusion 
equation, these higher-order contributions seem to reflect the changes needed 
for stable discretisations of equations with large amounts of advection, Mj., or 
dispersion, Uxxx- The extra complications of these nontrivial effects of odd 
derivatives appear necessary. However, here we restrict attention to proving 
consistency to an error quadratic in the odd coefficients and leave to further 
research the investigation of higher-order consistency. 

Theorem 3 The centre manifold model ^ constructed with the amplitude 
condition the internal boundary condition ^) and to errors 0{^^,e^^, 
forms a finite difference approximation to the PDE Ut = Au + eBu consis- 
tent to O (^d1^ + ed1^~^^ e^j, where A is any even operator and B is any odd 
operator. 

Proof: As in pO|), we expand the the centre manifold ansatz (|^ to errors 
u = Y.^^v' + eY.^'w\ and = X: 7 + 6 E , (32) 

fc=0 fc=0 fc=0 fc=0 

^One easily imagines other functions pk and that have all the requisite properties 
to ensure a consistent finite difference approximation g{u). For example, PkiC) = Pk{0 + 
asin(2n7r^) . However, as is often the case, if the method of construction of the centre 
manifold makes v'^ a polynomial of degree 2k, then the solution for w'^ is the one given 
in (|2^) in terms of pk and ■ 
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where and are some difference operators, and and are functions 
defined about the jth gridpoint, and they all implicitly refer to the jth el- 
ement (as before the superscript to v, w, g and h denotes an index in the 
series, whereas the superscript to 7 denotes exponentiation). After substi- 
tution into the pde, terms in e° determines and g'^ as in the previous 
section. Upon extracting from the pde the coefficients of terms linear in e 
and of various powers of 7 requires us to solve 

k 

Aw'' + Bv'' = Y^ [f-'-v' + g^'-^'w') , for < A; < £ . (33) 

Substitution of the expansion (^) into the amplitude condition (0) and the 
iBC's (^ and equating coefficients of £7^ leads to these internal boundary 
conditions for the w''{x): 

w'^ = for X = Xj,Xj±i . (34) 

Since B is an odd operator, we formally write it as the following infinite sum 
of odd powers of centred difference operators 

00 

B=Y: hmliJl^-' , (35) 

m=l 

for some coefficients 6^; for example (from [0, p65,e.g.]) 

-— = — — arcsmh(;rOa.) = —UxOx tI^xK H r^-xK — ■ ■ ■ 

and more generally, B could be any antisymmetric convolution operator for 
which the infinite sum ( P^ ) is reasonable. I prove that there exists solutions 
w'^(x), odd functions of x (about Xj), with 

f = hfxS^''-' , (36) 

so that the model ( p2| ) is consistent with the effect of the odd derivatives in 
eB to errors O (S^^^^ = O {d1^~^^. Since we already know v'^ and since w'^ 
appears to vary for different problems, the operators are determined by 
the solvability condition that all terms except Aw'' — g^w'' appearing in (|33|) 
combine to be in the range of the singular A — g^. 

• First, prove that the even part of Bv'' cancels with the even part of 
Y^'^=Q f^~^ and so is eliminated from (^3]). As a preliminary step 
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consider, using (|3lD, 

1 1 



{2k-iy.2 
1 



k—1 k—1 

n i^+i-m)- n (e-i-"^) 

m=— fe+l m=—k+l 



(2A; - 1)(2A; - 2) 
= (37) 

Then from (|^) and since pk is odd and is even, see Figure Bqk is 
odd and so the even part of Bv'^ is 



Bpkfid^'-'u, = J2b^fiJl"'-'pkfi6"'"\ by (35) 

m=l 
k 

= Yl bmfixSxPk-m+1 fi6^'''^Uj by d^) inductively 

m=l 

k 

= Y.bmqk-mf^S^''-'u, by(^. (38) 

m=l 

Whereas on the right-hand side of (|33D the even part of Y.r=o f''~^v^ is 

r=0 

j^q'bk-rfiS^'-'u,, (39) 



r=0 r=0 



r=0 



which exactly cancels with ( |57| ) from the even part of i3f on the left- 
hand side of (BSl). 



Second, since the even components of ( |5BD that involve the various v'^ 
cancel, and since A and g'' are even, then a particular solution w^{x) 
of (|33[) may be found that is odd. The conditions (|3^ ) then force 
these odd w'^ to be the unique solutions in the space of finite degree 
polynomials. 

Third, consider evaluating the hierarchy of equations ( P^D at the centre 
grid point of the element, x = Xj or equivalently ^ = 0, and simplify 
the various contributions. 
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— Now A is an even operator and an odd function, so Aw^ is an 
odd function, and thus Aw''{xj) = 0. 

— Since g''~^ is a difference operator in j it does not affect the am- 
phtude condition that w^{xj) = 0, and thus g''~^w^{xj) = 0. 

— I have already shown that the even parts of i3f ^ and Y,r=o f^~^v^ 
agree pointwise, so they certainly do at Xj, at which the odd parts 
must also trivially vanish together. 

Thus the choice (^6]) is the unique one to satisfy this solvability condi- 
tion for (^). 

Since the expansion ( |5^ ) then contains the exact differences up to bi_ifi6^^^^ 
and a£_i5^^~^, the errors in the finite truncation of the finite difference model 
will be O (d'^^^ from the even terms and O (^ed'^^~^, e^^ from the odd. 

Corollary 4 It follows immediately from the theorem that if the highest 
derivative in the operator A + eB is of order n, then the finite difference 
model i^dj) is consistent with Ut = Au + eBu to O (ji'^^^^^'^^ as h ^ to an 

error O {e^). 



6 Concluding remarks 

We have seen that the artificial internal boundary conditions together 
with the application of centre manifold theory generate finite difference mod- 
els that have remarkably good properties, at least for linear systems. These 
are explcitly shown for the example of the advection-diffusion equation 
where we saw not only consistency for small h, but also an appropriate up- 
wind model for large eh. Although the Theorem ^ only establishes consis- 
tency with the odd terms to O (e^) the advection-diffusion example of §0 
shows that higher-order consistency is possible. Further research is needed 
on the characteristics of higher orders in the odd derivatives. 



Also, further research, such as that for Burgers' equation in will 
explore the performance of this holistic approach to discretisations of non- 
linear systems in various spatial dimensions. Since centre manifold theory is 
designed to analyse nonlinear systems I expect reliable models to be derived. 

Throughout the analysis in this paper I have parametrised the centre 
manifold model in terms of the field at each of the grid points, Uj = u{xj, t). 
This was done for simplicity. Other choices are possible for the parameters 
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of the finite difference model, for example we could choose to use the element 
average 

^ii^) = T / u{x,t) dx. (40) 

h Jij-/i/2 

This choice would be appropriate to easily establish the conservation of total 
or not as the case may be. Computational experiments show that either 
of these choices of amplitude produce equivalent results for linear systems. 

Centre manifold theory is routinely applied to autonomous dynamical 
systems. However, the geometric viewpoint it establishes leads to a rational 
treatment of the projection of initial conditions onto the finite dimensional 
model and of the projection of a perturbing forcing [^, ^ ^]. 

Acknowledgement: this research was supported by a grant from the Aus- 
tralian Research Council. 



A Computer algebra for perturbed diffusion 

This computer algebra program is included to replace the recording of te- 
dious details of elementary algebra involved in solving the advection-diffusion 
PDE (ISl). 



For this problem the iterative algorithm of [IS] is implemented in RE- 
duce]^ Although there are many details in the program, the correctness of 
the results are only determined by driving to zero (line 40) the residuals of 
the governing PDE, evaluated on line 31, and the residuals of the iBC's, eval- 
uated on lines 32-35, to the error specified on line 29. The other details in 
the program only affect the rate of convergence to the ultimate answer. 

Lines 44-46 then rewrite the model in terms of the central difference 
operator 5 and mean /i. Lastly, lines 48-53 derive the equivalent differential 
equation for the finite difference model. 

Other problems in the same class can be and have been examined by 
amending line 31, the evaluation of the residual of the governing PDE, with 
other differential operators such as 

+eps* (a0*v+a2*df (v , x , 4) +b2*df (v , x , 3) ) 

The results of the analysis with such terms confirm the theorems in Sections ^ 
and H 

^At the time of writing, information about reduce was available from Anthony 
C. Hearn, RAND, Santa Monica, CA 90407-2138, USA. [mailto : reduceQrand . org| 
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1 % improve printing 

2 on div; off allfac; on revpri; factor gam; 
3 

4 % make function of xi=(x-x_j)/h 

5 depend xi,x; 

6 let df (xi,x)=>l/h; 
7 

8 % parametrise with evolving u(j) 

9 operator u; depend u,t; 

10 let df (u(~k) ,t)=>sub(j=k,g) ; 
11 

12 % solvability condition 

13 operator solg; linear solg; 

14 let { solg(xi-~p,xi)=>(l+(-l)~p)/(p+2)/(p+l) 

15 , solg(xi,xi)=>0, solgd ,xi)=>l }; 
16 

17 % solves v"=RHS s.t. v(0)=0 and v(+l)=v(-l) 

18 operator solv; linear solv; 

19 let { solv(xi"~p,xi) => 

20 ( xi-(p+2)-(l-(-l)-p)*xi/2 )/(p+l)/(p+2) 

21 , solv(xi,xi) => (xi"3-xi)/6 

22 , solv(l,xi) => (xi''2)/2 }; 
23 

24 °/o linear solution in jth interval 

25 v:=u(j); 

26 g:=0; 
27 

28 % iterative refinement to specified error 

29 let { gam"3=>0, eps~4=>0 }; 

30 repeat begin 

31 deq:=-df (v,t)-eps*df (v,x)+df (v,x,2) ; 

32 abc : =- (sub (xi=+l , v) -sub (xi=- 1 , v) ) /2 

33 +gam*(u(j+l)-u(j-l))/2; 

34 bbc : =- (sub (xi=l , v) -2*sub (xi=0 , v) +sub (xi=-l , v) ) 

35 +gam*(u(j+l)-2*u(j)+u(j-l)) ; 

36 gd:=bbc/h''2+solg(deq,xi) ; 

37 g:=g+gd; 

38 V : =v+h"2*solv(-deq+gd , xi) +xi*abc ; 

39 showtime; 

40 end until (deq=0)and(abc=0)and(bbc=0) ; 

41 deq:=sub(xi=0,v)-u(j) ; % confirms amplitude 
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42 

43 % get difference form of evolution 

44 gd:=( g where { mu"2=>l+del"2/4, u(j)=>l, 

45 u(j+~k)=>(l+mu*del+del~2/2)"k when k>0, 

46 u(j+~k)=>(l-mu*del+del~2/2)"(-k) when k<0 } ) ; 

47 % deduce the equivalent differential form, 

48 operator du; 

49 o:=8; 

50 let h~~p=>0 when p>o; 

51 ge:=( g where u(~k)=>(du(0)+f or n:=l:o+3 sum 

52 du(n)*h''n*(k-j)'"n/factorial(n)) )$ 

53 gee:=coeff (sub(gam=l,ge) ,h) ; 

54 

55 end; 
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